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1  Introduction 


Fast  Multipole  Method  (FMM)  belongs  to  a  class  of  algorithms  developed  for  the  rapid 
evaluation  of  Coulombic  interactions  in  large-scale  particle  ensembles.  In  two  dimen¬ 
sions,  the  existing  implementations  of  the  FMM  rely  on  Taylor  or  Laurent  (multipole) 
expansions  for  the  potential  field  (see  [2], [4]),  while  the  three-dimensional  ones  are  based 
on  spherical  harmonics  (see  [4]).  The  expansions,  as  well  as  the  corresponding  transla¬ 
tion  operators,  are  obtained  algebraically  from  explicit  formulas  relating  the  Newtonian 
potential  and  its  partial  derivatives  at  various  locations. 

In  this  paper,  we  construct  a  version  of  the  FMM  based  on  a  different  analytical 
apparatus.  Instead  of  multipole  and  Taylor  expansions,  we  use  specially  designed  bases, 
consisting  of  singular  functions  of  an  appropriately  chosen  operator.  The  expansions 
we  use  display  much  faster  convergence  than  the  previously  used  ones.  In  addition, 
we  introduce  an  intermediate  representation  consisting  of  complex  exponentials,  and 
diagonalizing  most  translation  operators.  When  these  two  techniques  are  combined,  the 
resulting  algorithm  is  about  five  times  faster  than  the  old  one  for  reasonably  uniform 
distributions,  and  about  three  times  faster  for  highly  non-uniform  ones. 

The  structure  of  this  paper  is  as  follows.  Section  2  introduces  the  mathematical 
preliminaries.  In  Section  3  we  describe  the  analytical  apparatus  to  be  used.  Section  4 
contains  a  detailed  description  of  the  algorithm,  together  with  its  complexity  analysis. 
Numerical  experiments  and  the  performance  of  the  scheme  are  discussed  in  Section  5. 


2  Mathematical  preliminaries 


2.1  Notation 

For  any  number  s  >  0  we  will  denote  by  Df  the  boundary  of  the  square  [-f ,  |]  x  [-|,  |], 
and  by  the  boundary  of  the  square  [-y,  f]  x  [-^,  f  ]  (see  Figure  1). 

We  will  call  fi’"  the  open  square  within  the  inner  square  £>*”  and  the  region 
outside  the  outer  square  If  s  =  1,  we  will  simply  write  Z)*”  for  Df,  Z)°“*  for 
n*"  for  and  for 

We  will  denote  by  the  set  consisting  of  the  rectangle  [-^,  x  [-=^,  y]  minus 

two  squares  [-f ,  x  [-f ,  and  [-f ,  f  ]  x  =^]  (see  Figure  2).  _ 

For  any  s  >  0  and  a  complex  number  Zq  =  xo  +  iyo,  we  denote  by  the  open _ 

square  (xq  -f,a:o  +  f )  x  {yo  -  2/o  +  §)•  Finally,  we  define  the  region  0°“*  by  Sie  formula  3 

Qout  ’  □ 

S^Zq  \  3s, Zq  *  Qj 

For  any  set  5  C  we  will  denote  its  closure  by  5.  _ 

In  agreement  with  standard  practice,  we  will  denote  by  P  the  Hilbert  space  of  all - 

complex  sequences  x  =  {x^},  such  that  |  ^nP  <  oo,  with  the  inner  product  defined  — 

1  I  Availsfe.iiijfe,y 

avail,  aai/ir"” 
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Figure  1:  The  domains 
and  their  boundaries. 


by  the  formula 


(1) 


(X,y)  =  Y^XuVn- 

n=l 

The  standard  basis  of  P  will  be  denoted  by  61,62,  - . .. 

Let  X  be  a  piecewise  smooth  curve  in  A  function  f  :  X  €  belongs  to  the 
vector  space  L^(X),  if  and  only  if, 

Wfhnx)  \f{x)\'^  dx^ "  <  00,  (2) 

where  the  integration  is  performed  with  respect  to  the  arclength  dx.  We  say  that  a 
function  f  :  X  €  belongs  to  the  space  L°°{X),  if  and  only  if, 

||/IU“(X)  ess  sup  \f{x)\  <  00.  (3) 

x&X 


2.2  Electrostatic  potentials  in  two  dimensions 

In  this  section  we  list  several  facts  from  mathematical  analysis,  which  will  be  used 
throughout  the  paper;  all  of  them  are  either  well-known,  or  follow  immediately  from 
well-known  results. 

A  unit  charge  located  at  the  point  Xo  €  generates  a  potential  aiid  a  field  given, 
respectively,  by  the  expressions 


*x.(x)  =  -log(||x-Xo||), 


£x.(x) 


X-Xq 

|x-Xo||2’ 


In  this  paper  we  will  work  with  analytic  functions  of  a  complex  variable,  making 
distinction  between  a  point  x=(x,z/)€E^  and  a  complex  number  z  =  x  +  iy. 

Since 


(4) 

(5) 


(x)  =  -Re  (log(2  -  2:0)),  (6) 

following  standard  practice,  we  will  refer  to  the  analytic  function  log(z)  as  the  poten¬ 
tial  due  to  a  charge.  To  describe  more  complicated  charge  distributions  we  will  need 
derivatives  of  log(2:),  and  we  will  also  refer  to  them  as  potentials. 

The  following  lemma  is  an  immediate  consequence  of  the  Cauchy-Riemann  equations. 


Lemma  2.1  If  f  :  G  G  is  analytic  and 


t.(z)  =  Re  (J(z)), 

then 

Wu  =  (uj,,  Uy)  =  {Re  (/),  -Im  (/')). 


(7) 

(8) 
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2.3  Hilbert-Schmidt  theory  for  integral  operators 


Let  X  and  Y  be  piecewise  smooth  curves  in  We  will  be  working  with  Lebesgue  spaces 
L^(X)  and  L^(y)  of  functions  square  integrable  with  respect  to  arclength  measures  dx 
and  dy.  The  product  space  L^(X  x  Y)  consists  of  functions  k  :  X  x  Y  ^  such  that 


r 

ll^||L2(;rxy)  =  JJ  dx  dy 


X-kY 


<  oo. 


(9) 


We  start  with  a  well-known  lemma  (see,  for  example,  [7],  sec.VI.6). 
Lemma  2.2  IfkE  L^{X  x  Y),  then  the  expression 

^/(®)  =  JyKx,y)f{y)dy 

defines  a  continuous  operator  A  :  L^{Y)  — ^  L^{X),  and 


(10) 


ll^ll  <  ||^||L2(jfxy)  •  (11) 

Remark  2.3  The  integral  operator  induced  by  a  kernel  k  €  L'^{X  x  Y)  via  formula  (10) 
is  usually  referred  to  as  a  Hilbert-Schmidt  operator  (see,  for  example,  [7]). 

The  following,  is  an  immediate  consequence  of  Theorem  VI.  17.  in  [7] 

Theorem  2.4  IfkE  L'^{X  X  Y),  then  there  exist  two  orthonormal  systems  of  functions 
in  nnd  {V’n}  in  L^{X),  and  a  sequence  {s„},  n=l,2,. . of  non-negative  real 

numbers  such  that 


1-  X)  (12) 

71=1 

OO  _ 

2-  k(^X,  y)  —  ^  ^  ^n(®)  ^n{,y^  5  (1^) 

n=l 

in  L^i^X  X  V)  sense.  Moreover,  the  sequence  {5,1}  is  uniquely  determined  by 
k  €  L^{X  X  Y). 


Theorem  2.5  (Canonical  form  for  Hilbert-Schmidt  operators) 

Let  A  :  L^{Y)  — >  L^{X)  be  the  integral  operator  induced  by  a  kernel  k  G  L'^{X  x  Y)  via 
expression  (10).  Then,  for  any  f  €  L^(Y), 

OO 

^/  =  S  (  /,  <f>n)  ,  (14) 

n=l 

where  the  functions  ipn,  <i>n  nnd  numbers  Sn  are  provided  by  Theorem  2.4-  Equivalenty, 


for  all  n  =  1,2,.. . . 


(15) 
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Remark  2.6  A  finite  dimensional  version  of  Theorem  2.4  is  known  in  numerical  analysis 
as  the  Singular  Value  Decomposition  (SVD);  the  coefficients  are  referred  to  as  the 
singular  values  of  the  operator  A,  and  the  functions  and  <t>n  are  called  left  and  right 
singular  vectors,  respectively. 

To  restate  Theorem  2.4  in  these  terms,  we  define  operators  U  :  P  and 

V  :  P  L'^iY)  by  specifying  them  on  the  elements  of  the  standard  basis  {e^}  in  P  (see 
Sec.  2.1),  via  the  formulae 

U Cfi  —  Ipn  ) 

VCn  =  <j>n  ■  (17) 

Now,  (14)  can  be  rewritten  in  the  form 

A  =  UDV\  (18) 

where  D  :  P  P  is  a,  diagonal  operator  Avith  the  coeflBcients  s„  on  the  diagonal,  so  that 
for  all  n  =  1,2,. . . , 

D^n  ~ 

As  in  the  finite  dimensional  case,  U  and  V  are  isometries. 

Remark  2.7  Given  an  operator  A  :  L'^{Y)  L^{X)  defined  by  the  formula  (10),  the 

operator  A^  :  L^(^X)  — >  L^(Y^  is  referred  to  cis  the  transpose  of  A  if  and  only  if  for  all 
/  €  L^iY)  and  g  €  L^{X) 

=  {A^gJ)-  (20) 

Similarly,  the  operator  A*  :  LP{X^  L^(Y)  is  referred  to  as  the  adjoint  of  A  if  and  only 
if  for  all  /  €  L^Y)  and  g  G  L^{X) 


{Af,g)  =  {f,A*g).  (21) 

The  following  well-known  lemma  gives  explicit  expressions  for  the  operators  A^  and  A* 
in  terms  of  the  kernel  of  A. 

Lemma  2.8  If  an  operator  A  :  P{Y)  -4  LP{X)  is  defined  by  (10),  then  the  operators 
A'^ , A*  :  L^{X)  —*  L^{Y)  are  defined,  respectively,  by  the  formulae 


^^/(y)  =  y)  fix)  dx, 

(22) 

^*fiy)=  f  k{x,y)  f{x)dx. 

•f  JC 

(23) 

The  following  lemma  follows  immediately  from  Lemma  2.8. 
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Lemma  2-9  Let  and  <l)n  be  the  left  and  right  singular  functions  the  operator  A  de- 

fined  by  (10),  and  $n  its  singular  values.  For  each  n=l,2,.. 
tpn  €  L‘^{Y),  (fy  €  L^{X)  via  the  formulae 

. ,  we  define  the  functions 

V’n  =  <l>n, 

(24) 

(26) 

Then  and  are,  respectively,  the  left  and  right  singular  functions  ofA^.  Moreover, 

(26) 

for  all  n=l,2, -  Similarly,  the  functions  ip*  €  L^{Y),  ^* 

n=l,2,. via  the  formulae 

€  L^{X)  defined  for  each 

V’n  =  <^n, 

(27) 

<l>n  = 

(28) 

are  left  and  right  singular  functions  of  A*,  respectively.  Furthermore, 

(29) 

for  all  n=l,2, _ 

3  Analytical  apparatus 

3.1  Efficient  representation  of  potentials 

We  will  define  the  operators  Cout  '■  L^(Z)’”)  and 

Cin  :  L2(D°“<)  ^  by  the  formulae 

c.„./(*2)  =  L  <*^(0. 

JDxn  ^  —  Z2 

(30) 

ds(a 

JDotit  Zi  —  ^ 

(31) 

where  the  integration  is  performed  with  respect  to  the  arclength  ds  (see  Section  2.1  for 
the  definitions  of  D’”  and  In  other  words,  the  kernel  k  :  x  Z)*"  — C  of  the 

operator  Cout  is  given  by  the  expression 


k{z2,zi)- — - — .  (32) 

Z1-Z2 

Clearly,  the  operators  Cout  and  Cin  satisfy  the  conditions  of  Theorem  2.4.  Moreover, 
Cin  =  (Couty  ■  Thus,  a  combination  of  Theorem  2.4  and  Lemma  2.9  leads  to  the  following 
result. 
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Theorem  3.1  There  exist  orthonormal  systems  and 

in  and  non-negative  real  numbers  s„  such  that 


=  V’r,  (33) 

€  =  V’rS  (34) 

^in^n  ~  '5nV’n”>  (35) 

C^4'T  =  (36) 

for  all  n=l,2, - Moreover, 

1  °° 

k{z2,  2i)  =  - - -  ='£sn  C(^l)  €“'(^2).  (37) 

“  ^2  „=1 


The  following  theorem  provides  estimates  for  the  magnitude  of  the  singular  values 
and  singular  function  of  the  operator  Cout  ’  T^(I>’”)  -»•  T^(i)®“‘).  Its  proof  is  somewhat 
involved,  and  can  be  found  in  Appendix  A. 

Theorem  3.2  There  exist  constants  0  <  <  9  <  1  and  c>  0,  ci  >  0  such  that 


CiQi 

< 

(38) 

Sn 

< 

cq  , 

(39) 

IICIU 

< 

cn^ 

(40) 

< 

cn^ 

(41) 

for  all  n=l,2,. . . . 

Remark  3.3  Our  numericcd  experiments  show  that  the  maximum  values  of  the  left 
singular  functions  are  uniformly  bounded,  while  left  singular  functions  grow  as 
logn.  Furthermore,  the  coefficient  q  in  (39)  is  less  then  0.37  (see  Table  1).  However,  the 
crude  estimates  (38-41)  are  sufficient  for  the  purposes  of  this  paper. 

The  following  theorem  states  that  the  left  singular  functions  of  the  operator 
Cout  are  restrictions  to  of  functions  analytic  on  0°“*,  which  are  bounded  at  infinity. 
Similarly,  the  left  singular  functions  of  the  operator  are  restrictions  to  D®”  of 
functions  analytic  on  fi®'"  (see  Section  2.1  for  definitions  of  T>®”,  fi®",  and  fi°®®‘). 

Theorem  3.4  Under  the  assumptions  of  Theorem  3.1,  for  each  n=l,2,...,  there  exist 


complex  analytic  functions  ; 

C  and  ^j;®  :  fi®'®®  C  such  that 

1. 

1 

D*»  =  V’n”  > 

(42) 

2. 

1  pout  —  tpn  1 

(43) 

3. 

lim  •iT‘{z)  =  0. 

2:  -+00 

(44) 
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Table  1:  Singular  values  s„  of  the  operators  Cin  and  Cout 


n 

Sfi 

^n+1 / 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

0.2359412633095143E+01 
0.8200026261260834E+00 
0.2969654528722792E+00 
0.1117065933565274E+00 
0.3899373405395550E-01 
0.1419762842492349E-01 
0.5167619274037499E-02 
0.1882316460599726E-02 
0.6806476236529151E-03 
0.2474684420248163E-03 
0.9002598391003066E-04 
0.3279166684309465E-04 
0.1189492698766364E-04 
0.4326403854315149E-05 
0.1574295791427948E-05 
0.5732869472275525E-06 
0.2081873704334665E-06 
0.75738741251441 14E-07 
0.27562317941 17075E-07 
0.1003533522800739E-07 
0.3646364216302235E-08 
0.1326702566316864E-08 
0.4828202607541905E-09 
0.1757742306275760E-09 
0.6388965188123934E-10 
0.2324723901442769E-10 
0.8460347483779406E-11 
0.3079818643909223E-11 
0.1119689433835681E-11 
0.4074314198593252E-12 
0. 1482764897026021 E-12 
0.5397405796606410E-13 
0.1962573461216843E-13 
0.7141555758338406E-14 
0.2599021334444648E-14 
0.9460284166056962E-15 
0.3440293867375512E-15 
0.1251897831151196E-15 

0.3475452384309655E+00 
0.3621518314828151E+00 
0.3761602310170769E+00 
0.3490728065576352E+00 
0.3641002527554371E+00 
0.3639776390376512E+00 
0.3642521557376765E+00 
0.3616010580049082E+00 
0.3635779123075419E+00 
0.3637877346033592E+00 
0.3642466921 090881E+00 
0.3627423712426653E+00 
0.3637184035515401E+00 
0.3638809145978704E+00 
0.3641545320448064E+00 
0.3631468873315051E+00 
0.363800844853 1035E+00 
0.3639130712466957E+00 
0.3640962000883560E+00 
0.363352507261 1107E+00 
0.3638425806137019E+00 
0.3639250220903513E+00 
0.3640572795205558E+00 
0.3634756451678426E+00 
0.3638654825924015E+00 
0.3639291306 1 75649E+00 
0.3640298049003309E+00 
0.3635569373703303E+00 
0.3638789538842048E+00 
0.3639299339108355E+00 
0.3640095478003274E+00 
0.3636142130448670E+00 
0.3638873091614348E+00 
0.3639292924948542E+00 
0.3639940942646556E+00 
0.3636565040740656E+00 
0.3638927020226409E+00 
0.3639280307828200E+00 
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Proof.  It  follows  from  (38)  that  >  0,  for  all  n=l,2,. . . .  Thus,  (15)  can  be  written  as 

=  s;’  c„,c‘(22)  =  -s;'  / .  7^  *({),  (45) 

Jl,,n  ^  —  Z2 

for  any  Z2  €  Therefore,  the  formula 

=  (46) 

extends  to  a  function  continuous  on  and  analytic  in 

In  a  similar  manner,  we  extend  xj>^  to  an  analytic  function  on  the  domain  fi®",  con¬ 
tinuous  to  the  boundary  D*”,  by  setting 

K"(^i)  =  Cj.C(^.)  =  ^  &(0-  (4T) 

□ 

Given  a  charge  anywhere  in  the  region  O’”,  the  following  two  theorems  allow  us  to 
represent  its  potential  inside  the  region  0°“*  by  a  linear  combination  of  left  singular 
functions  of  the  operator  Covt-  Similarly,  the  potential  of  a  charge  in  the  region  0°”* 
can  be  expressed  inside  the  region  O’”  as  a  linear  combination  of  left  singular  functions 
of  the  operator  C,„. 

Theorem  3.5  For  any  zi  e  P’”,  Z2  €  D®”*, 

1  °° 

=  E^nC(^l)V^r(^2).  (48) 

n=l 

Moreover,  there  exists  a  constant  ci  >  0,  such  that  for  any  zi  €  Z2  €  and 
integer  iV  >  0 

I  -  V  -  £  ’‘?‘'(^2)  \<c,N^  (49) 

“  •^2 

Proof.  Let  us  define  Sn,  N=1,2,.  . . ,  by  the  formula 

Sat  =  f;  (50) 

n=l 

It  follows  immediately  from  Theorem  3.2,  that  for  there  exists  a  constant  ci  >  0,  such 
that  for  any  integers  M  >  N  >  0, 

M  M 

=  II  E  ^.«”«“'IU<  E  «»ll«”ll»liv>riu 

n=iV-fl  n—N+1 

<  f;  c^n\^  <ciN^q^.  (51) 

n=iV+l 
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Since  CiN'^q^  0  as  N  oo,  the  sequence  {5jv}  of  continuous  functions  En=l 
converges  in  the  maximum  norm  to  a  limit,  which  we  denote  by  S.  Due  to  Theorem  2.4, 
{S'iv}  converges  to  the  function  in  x  jD*”),  so  we  have 

S{z2,zi)  =  — - — ,  (52) 

Zi  —  Z2 

for  all  zi  G  D*'^,  Z2  G  Thus,  from  (51)  and  (52)  we  obtain 

II  — - E  «nC  V-riU  =  115  -  5k|U  =  II  lim  Sm  -  5jv|U  <  0i  (53) 

Zl  —  Z2  M-^oo 

□ 


Theorem  3.6  For  any  zi  G  D*",  Z2  G  0°“*, 

1  °° 

—  =E"n^r(^l)^r(^2). 


Zl  -  Z2 


(54) 


n=l 


Furthermore,  there  exists  a  constant  Ci  >  0,  such  that  for  any  Zj  G  D*”,  Z2  G  and 

integer  N  >  0, 

I  -  E  I  <  c.  jv^ (55) 

-  ^2  ri=l 

Proof.  Due  to  a  combination  of  the  maximum  modulus  principle  for  complex  analytic 
functions  and  (49),  we  have 


1 


N 


E  I  = 


n=l 


max  max  I 

Zl  —  Z2 

I  7^  “  (^l)  ^n“*(^2)  I  = 

zieD'^  Z2eD^^*  Zl  —  Z2 
,  1 


max  max 


-  E  (^2)  I  <  CxiV^A 


(56) 

□ 


zieD'”  z2eD°'^t  Z1  —  Z2 

which  proves  (55),  and  (54)  follows  from  (55)  immediately. 

Let  t)e  the  fimctions  provided  by  Theorem  3.4.  For  any  real  num- 

r  5  >  0  and  point  zq  G  we  define  i 
^Ts,zo  '  ^°s%  n=l,2,. . . ,  by  the  formulae 

1  -  ^0 


ber  5  >  0  and  point  zq  G  we  define  analytic  functions  ^^n,s,zo  ■  ^l?zo  ® 


^^.....(^i)  =  ;;;;5  «r( 


), 


(57) 
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(58) 


for  all  zi  €  and  all  Z2  E  Moreover,  suppose  that  and  are  the 

functions  provided  by  Theorem  3.1.  For  any  real  number  5  >  0  and  point  zq  €  D®",  we 

define  functions  0^"  ^  C  and  ^  <D,  n=l,2,. . . ,  by  the  formulae 

^  «*(^^).  (59) 

KLM  =  v-r‘(^^).  (60) 

for  all  zi  €  and  all  Z2  €  B^X,  (see  Section  2.1  for  definitions  of  Z)®’\ 

Qin  Qout  \  *  *  ’  ’ 

^5,20 5^0 /* 

The  following  theorem  is  an  immediate  consequence  of  Theorem  3.6  and  Theorem  3.1. 

Theorem  3.7  Suppose  that  s  >  0  is  a  real  number  and  zq  E  D*®®.  Then 

1.  The  functions  {V’n^s.zol  form  an  orthonormal  system  in  T^(Z)®"^^).  The  functions 

{^n!s,2o}  form  an  orthonormal  system  in  L^(BfX)- 

O  ®out  _  ^out  „  Y  g 

^71,5,20  I  DOW*  ^71,5,20^ 

3.  For  any  zx  €  Z2  E 

=  E  ^n»t,..{z.)  (61) 

~  ^2  n=l 


Furthermore,  there  exist  constants  c>0,  0<9<1  such  that  for  any  s  >  0,  zi  £ 
Z2  E  ,  and  integer  N  >  0, 


1 


N 


r.N 


^1-^2 

The  following  two  theorems  axe  immediate  consequences  of  Theorem  3.7. 


(62) 


Theorem  3,8  Suppose  that 

m 

«(-)  =  E-rT  (63) 

t=l 

is  the  potential  due  to  a  set  of  m  charges  of  strengths  91,92, •••,9m  located  at 
points  Zx,  Z2., . . . ,  Zm  inside  the  square  Suppose  further  that  the  functions 

Ks,zo  •  ^Zo  KZo  •  Kl  ®  h  the  formulae  (57),  (58),  and 

are  defined  by  (37).  Then  for  any  z  €  0^* , 


*(z)=i:<.„«~‘^(2), 


n=l 


(64) 
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with 


On  —  •Sn  ^n!s,zo 

i=l 

for  all  n  =  1,2, -  Furthermore,  there  exist  constants  c  >  0  and  0  <  g  <  1,  such  that 

for  any  real  number  s  >  0,  integer  N  >  0  and  z  G  , 

i*(z)  -  f;  i,ii.  (66) 

Theorem  3.9  Suppose  that  $  given  by  the  formula  (63)  is  the  potential  due  to  a  set  of 
m  charges  of  strengths  9i,  92,  •  •  • ,  9m  located  at  points  zi,Z2,...,Zm  inside  the  region  . 
Suppose  further  that  the  functions  ^°n%o  '  ®  defined  by 

the  formulae  (57),  (58),  and  s„  are  defined  by  (37).  Then  for  any  z  € 

m  =  -£a.¥Z^{z),  (67) 

n=l 

with 

m 

On  =  -5ni;9«’^Wo(^»)>  (68) 

i=l 

for  all  n  =  1,2, -  Furthermore,  there  exist  constants  c  >  0  and  0  <  q  <1,  such  that 

for  any  real  number  3  >  0  integer  N  >0  and  z  € 

|*(^)  -  (69) 

n=l  ^ 

3.2  IVanslation  Operators  and  Error  Bounds 

The  following  five  theorems  allow  us  to  translate  expansions  of  the  forms  (64),  (67)  from 
one  center  zq  to  another,  and  to  convert  expansions  of  the  form  (64)  into  expansions  of 
the  form  (67).  We  only  provide  proofs  of  Theorems  3.10,  3.11,  3.16  below;  the  proofs  of 
Theorems  3.12,  3.14  are  virtually  identical  to  the  proof  of  Theorem  3.10,  while  the  proofs 
of  Theorems  3.13,  3.15  are  identical  to  the  proof  of  Theorem  3.11.  Thus,  the  proofs  of 
Theorems  3.12,  3.13,  3.14,  3.15  axe  omitted. 

For  a  real  number  s  >  0  and  a  point  zq  G  we  define  coefficients  p°rt],^s,zo  ^^6 
P^n,k,s,zo  ’  respectively,  by  the  formulae 

(70) 

(71) 
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If 

II 

’([fOUt 

_ 

/. 

‘'■^23,0 

=  1,2,.... 

Theorem  3.10  Suppose  that  #  defined  by  the  formula  (63)  is  the  potential  due  to  a 
set  of  m  charges  of  strengths  91,92,  •••,9m  located  in  the  square  and  are 

functions  defined  by  (58),  with  n  =  1,2, - Suppose  further  that 

=  (72) 

71=1 

00 

$(z)  =  •£  (73) 

n=l 

are  the  expansions  provided  by  Theorem  3.8  valid  in  and  respectively.  Then 

for  n  =  1,2,... 

=  (74) 

*=1 

with  the  coefficients  p^k,s,zo  defined  by  (70). 


Proof.  Since  the  functions  form  an  orthonormal  system  in 

r  _ 

=  L.  E‘‘'‘S.,o(^)«.,o(^)*(^) 

'^^25,0  A:=l 

00  _  00 

=  L.  =  EpS, 

''^is.O  fc=l  '  - 

for  n  =  1,2,. . . . 


A:=l 


(75) 

(76) 
□ 


Theorem  3.11  Suppose  that  under  the  conditions  of  Theorem  3.10,  N  >  0  is  an  integer 
and  coefficients  b((  (n  =  1,2,.  ..,N),  are  defined  by  the  formula 

bn  (77) 

k=l 

K?fc,5,2o  defined  by  (70).  Then  there  exist  constants  c>  0  and  0  <  9  <  1,  such  that 
for  any  s  >  0,  N  >  0  and  z  € 

N  -  m 

l»(^)  -  E  <  ;  JV*,"  J:  |«|.  (78) 

n=l  ^  i=l 

Proof,  Let  0  <  9  <  1  be  the  constant  provided  by  Theorem  3.8.  Clearly,  it  is  sufficient 
to  show  that 


N 


N 


E».®S.W-E*^'«:5,oWI  <  -NYj:ki\, 


71=1 


71=1 


(79) 


1=1 
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for  some  constant  c>  0.  Now,  by  a  combination  of  (39),  (40),  (41)  and  (65), 


N 


N 


I E  ’’••KM  -  E  »S,o(^)l  =  I  E(‘«  - 

n=l  n=l  n=l 

oo  N  oo  m 

=  lE^rioW  E  ^  E^I  Tr  E  KlsJ^k-7=^\qi 

n=l  *=^+1  n=l  V-®  k=N+l  t=l 


_  7n  N  <x> 

^  7^“El«il  E  E 

*  2=1  n=l  fc=N+l 


(80) 


for  some  constants  c\  >  0,  C2  >  0.  Due  to  Schwarz’s  inequality  and  (41),  there  exists 
C3  >  0,  such  that  for  all  positive  integers  k  and  n 


\P^k,s,zo  \  ^  ^  C3/;. 

Now  (79)  follows  from  (80)  and  (81). 


(81) 

□ 


Theorem  3.12  Suppose  that  #  defined  by  the  formula  (63)  is  the  potential  due  to  a 
set  of  m  charges  of  strengths  qi,q2,. . .  ,qm  located  in  the  square  0°^^,  and  '^^n^s,zo 
functions  defined  by  (57),  with  n  =  1,2, - Suppose  further  that 

00 

*{z)  =  y;a.«“.,„W  (82) 

n=l 

00 

»(^)  =  E^»«t,..(^).  (83) 

n=l 

are  the  expansions  provided  by  Theorem  3.9  valid  in  O^s.o  respectively.  Then 

for  n  =  1,2,.. .  formula 

OO 

^n  =  EK"wo«fc.  (84) 

fc=l 

with  the  coefficients  defined  by  (71). 

Theorem  3.13  Suppose  that  under  the  conditions  of  Theorem  3.12,  N  >  0  is  an  integer 
and  coefficients  bff  (n  =  1,2,. .  .,N),  are  defined  by  the  formula 

C  =  (85) 

fc=l 

with  p^)fk,s,zo  defined  by  (71).  Then  there  exist  constants  c  >  0  and  0  <  q  <  1,  such  that 
for  any  s  >  0,  integer  N  >  0  and  z  e  Offs) 


Suppose  that  s  >  0  is  a  real  number  and  zq  €  For  k  =  1,2,. . . ,  and  n  =  1,2,. . . , 
we  define  coefficients  qn,k,s,zo  by  the  formula 

MO-  (87) 

Theorem  3.14  Suppose  that  $  defined  by  the  formula  (63)  is  the  potential  due  to  a  set 
ofm  charges  of  strengths  qi,q2,  ■  •  ■  ,qm  located  in  the  square  and 
functions  defined  by  (58)  and  (57),  respectively.  Suppose  further  that 

OO 

*(^)  =  £  (88) 

n=l 


*(^)= (89) 

n=l 

are  the  expansions  provided  by  Theorem  3.8  and  Theorem  3.9  valid  in  and  , 
respectively.  Then  for  all  n  =  lj2j _ the  formula 


OO 

9n,fc,5,2o  Ok, 

k=l 


with  the  coefficients  qn,k,s,zo  defined  by  (87). 


(90) 


Theorem  3.15  Suppose  that  under  the  conditions  of  Theorem  3.10,  N  >  D  is  an  integer 
and  coefficients  b^  (n  =  1,2,. . .  ,N),  are  defined  by  the  formula 


N 

~  ^^Qn,k,s,zo  Ok, 
k=l 


(91) 


with  qn,k,s,zo  defined  by  (87).  Then  there  exist  constants  c  >  0  and  0  <  ^  <  1,  such  that 
for  any  s  >  0,  integer  N  >  0  and  z  €  0,^*, 


N 


71=1 


(92) 


The  following  two  theorems  list  certain  properties  of  the  coefficients 
j  .  .1  1  ,  ,  ^  t'n,k,s,zoi  i'n,k,s,zo 

and  qn,k,s,zo’  Ihey  are  quite  similar,  and  we  only  prove  the  first  one. 


Theorem  3.16  Suppose  that  s  >  0  is  a  real  number,  zq  G  fi*”,  and  the  coefficients 
P^k,s,zo>  Pn,k,s,zo  defined  by  formulae  (70)  and  (71).  Then 


r'  n,/:,5,«o 

Pn,k,l,^^ 

(93) 

rJn 

Jrn,kfSyZQ 

_ 

(94) 

out 

“n,/!,S,2o 

—  D*" 

irk^n.,s^ZQ  7 

(95) 
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with  any  k  =  1,2,. ..,  n  =  1,2, _ 


Proof.  Due  to  (58)  and  (70) 


=  L.  K'H.AO  MO 

X-Zq.  1 


5  ' 

Zo.  1 


4 


put 


=  Pk,n,l,^^ 

which  proves  (93).  The  proof  of  (94)  is  identical. 

Turning  our  attention  to  (95),  we  observe  that  due  to  Theorem  3.1, 


Pk,n,l,^^ 

and  (95)  follows  directly  from  (93),  (94)  and  (97). 


(96) 


(97) 

□ 


Theorem  3.17  Suppose  that  s  >  0  is  a  real  number,  zq  is  an  arbitrary  point  in 
and  the  coefficients  qn,k,s,zo  <^'cc  defined  by  formula  (87).  Then 

qn,k,s,zo  =  9n,fc,l,^>  (98) 

for  all  k  =  1,2,...,  n  =  1,2, _ 


3.3  Diagonal  Form  of  Translation  Operators 

In  this  section  we  construct  a  representation  of  potentials  (63)  in  which  the  translation 
operators  are  diagonal.  We  start  with  an  obvious  lemma. 


Lemma  3.18  If  z  and  zq  are  complex  numbers  such  that  Re  (z  —  zq)  >  0,  then 


1 


Z  —  Zo 


^-T(z-zo) 


(99) 
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Table  2:  Quadrature  nodes  and  weights  for  N  =  8. 


k 

Xk 

Wk 

1 

0.9743818326893713-01 

0.2469944279808820+00 

2 

0.4917660712345750-h00 

0.5328649552527160+00 

3 

0.1147630094197554-f01 

0.7724853603957714+00 

4 

0.2028186101391804-F01 

0.9870187246041202+00 

5 

0.3121662089244688+01 

0.1200121694721202+01 

6 

0.4425614160158614+01 

0.1405423586612871+01 

7 

0.5925421085324173+01 

0.1592796104486006+01 

8 

0.7601289009353128+01 

0.1765198353876860+01 

One  of  principal  numerical  tools  of  this  paper  is  finite  quadratures  for  the  integrals 
of  the  form  (99),  approximating  them  by  expressions  of  the  form 

TOO  ^ 

J  (100) 

with  Wj,Xj  chosen  to  minimize  the  error  of  the  approximation.  It  turns  out  that  the 
classical  Laguerre  quadrature  requires  56  nodes  to  obtain  a  full  double  precision  (15-digit) 
approximation  to  (99)  for  all  zq  e  fi’i"  and  z  e  Bi  (see  Section  2.1  for  the  definition  of 
Bi);  it  requires  28  nodes  for  single  precision  (7-digit)  approximation  and  14  nodes  for 
the  3-digit  approximation.  In  this  paper,  we  use  quadratures  for  integrals  of  the  form 
(99)  designed  in  [9].  The  nodes  and  weights  for  the  quadratures  are  listed  in  Tables  2-5 
and  the  following  lemma  (proved  in  [9])  describes  the  performance  of  these  quadratures 
when  zo  G  z  E  Bi. 

Lemma  3.19  1.  If  the  nodes  X\,X2, .,xs  and  the  weights  'Wi,W2, . . . ,  tns  a,re  those  given 
in  Table  2,  then 

\— - <  10-®  (101) 

^0  k=l 

for  all  Zq  G  -2  €  jBi. 

2.  If  the  nodes  Xi,  X2, . . . ,  xjo  and  the  weights  Wi,  W2, . . . ,  tuio  are  those  given  in  Table  3, 
then 

1 

I - (102) 

^0  k=l 

for  all  Zq  €  fij”,  z  e  Bi. 

3.  If  the  nodes  Xi,X2,...,  xie  and  the  weights  Wi,W2, . . . ,  Wie  are  those  given  in  Table  4, 
then 

1 

I - <  10“’’  (103) 

^-^0  k=l 
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Table  3:  Quadrature  nodes  and  weights  for  N  =  10. 


_ ^ _ m _ 

0.79400973700479499E-01  0.20213268247442060E+00 
0.4059967502704461 7E+00  0.4452920131 0708538E+00 
0.95860548270566906E+00  0.65492570079022383E+00 
0.17076338623411169E+01  0.83991908942837771E+00 
0.26342522431201576E+01  0.10125227869573986E+01 
0.37330678114549477E+01  0.11856981580215332E+01 
0.50056635563091912E+01  0.13587490932348730E+01 
0.64476147019688304E+01  0.15237759923040743E+01 
0.80499560865687445E+01  0.16815303253729583E+01 
0.98062704155363723E+01  0.18393633495134454E+01 


Table  4:  Quadrature  nodes  and  weights  for  N  =  16. 

_ _ ■m _ 

0.50985651676060708E-01  0.13042997567403943E+00 
0.26521044170119003E+00  0.29617079732513146E+00 
0.63877523197814055E+00  0.44851876604376186E+00 
0.11575210696261956E+01  0.58677848054269013E+00 
0.18084194569861214E+01  0.71331490632911526E+00 
0.25812535628463790E+01  0.83115583049433043E+00 
0.34688375582342864E+01  0.94332332670457377E+00 
0.44670816002182331E+01  0.10530465654290459E+01 
0.55752269484049922E+01  0.11635980151835418E+01 
0.67950948336753639E+01  0.12764875588622437E+01 
0.81285324158592802E+01  0.1390341 7411114886E+01 
0.95753440067592314E+01  0.15029785005709016E+01 
0.11133967001554676E+02  0.16142090491898996E+01 
0.12804196003878283E+02  0.17271419059134570E+01 
0.14590645231278096E+02  0.18490423387442165E+01 
0.16505707680646142E+02  0.19928836723099878E+01 
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Table  5:  Quadrature  nodes  and  weights  for  N  =  33. 


k 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 


_ _ 

0.25222297727341369E-01 

0. 13247759128084677E+00 

0.32380669248023823E+00 

0.59671535823233030E+00 

0.94802704428024797E+00 

0.13742659894203061E+01 

0.18719475451288249E+01 

0.24377774897497581E+01 

0.30687681826830449E+01 

0.37622927481697938E+01 

0.45161007854489040E+01 

0.53283153101708052E+01 

0.61974248981676391E+01 

0.71222793332370721E+01 

0.81020919953769825E+01 

0.91364474341994351E+01 

0.10225307113296538E+02 

0.11368999161839711E+02 

0.12568170717361438E+02 

0.13823683086431786E+02 

0.15136452987314821E+02 

0.16507283044647311E+02 

0.17936750445208845E+02 

0.19425207754170005E+02 

0.20972906955216222E+02 

0.22580224099961821E+02 

0.24247959886492881E+02 

0.25977718733384464E+02 

0.27772422510974588E+02 

0.296371 18447457170E+02 

0.31580471971364523E+02 

0.33617895925272799E+02 

0.35778138778825402E+02 


_ Wk _ 

0.64680274571 235485E-01 
0.14962720607997425E+00 
0.23260602282627724E+00 
0.31267567589062277E+00 
0.38936352755688915E+00 
0.46253028974011126E+00 
0.53228190405863132E+00 
0.59887920022585347E+00 
0.66266415869490855E+00 
0.72400998959978491E+00 
0.78329340463991183E+00 
0.84088382086631104E+00 
0.89714360915324316E+00 
0.95243401551568258E+00 
0.10071219382388539E+01 
0.10615826335963539E+01 
0.11161921501942585E+01 
0.11713015487646939E+01 
0.12271878679241548E+01 
0.12839926657455333E+01 
0.13416829044598757E+01 
0. 14000734153500639E+01 
0.14589180106367176E+01 
0. 15180347903942609E+01 
0.15774212769044868E+01 
0.16373402708097615E+01 
0.16983905310960812E+01 
0.17616017849285610E+01 
0. 18286259025840425E+01 
0.19021749030833896E+01 
0.19870826536338642E+01 
0.20930613125138042E+01 
0.22428277080254514E+01 
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for  all  zq  €  z  E  Bi. 

4-  If  the  nodes  Xi,  X2, . . . ,  X33  and  the  weights  wi,W2, . . . ,  W33  are  those  given  in  Table  5, 
then 

1  33 

I - 10-^5  (104) 

^-^0  k=i 


for  all  Zq  €  fii";  z  £  Bi. 


3.4  Informal  Description 

Obviously,  approximating  the  integral  (99)  by  a  finite  quadrature  formula  gives  rise  to 
an  approximation  of  the  function  by  a  finite  linear  combination  of  exponentials  in 
each  of  the  regions  B\.  Indeed,  suppose  that  iVo  is  a  natural  number,  and  positive 
real  numbers  Xi,  X2, . . . ,  xjvo »  u^i,  1^2,  •  •  • ,  are  such  that 

TOO  -^0 

I  /  dx-Y,y^k  I  <  £,  (105) 

k-i 


for  all  Zq  €  ^  ^  ^  ^  sufficiently  small  positive  number.  Clearly,  (105)  can  be 

rewritten  in  the  form 


1 


Z  —  Zq 


No 

k=l 


(106) 


and,  given  an  arbitrary  point  w  €  Bi,  (106)  can  be  rewritten  in  the  form 


1 


Z  —  Zq 


No 

k=i 


(107) 


In  other  words,  the  potential  of  a  charge  at  the  point  zq  has  been  approximated  by  a 
linear  combination  of  exponentials.  For  a  potential  #  given  by  the  formula 


m 


»(^)  =  E 

t=l 


(108) 


due  to  a  set  of  m  charges  of  strengths  91, 92,  •  •  - ,  9m  located  at  the  points  zi,  22,  •  •  ■ ,  ^m  in 
the  square  flj",  the  expression  (107)  assumes  the  form 


m  No 

»(^)  =  E-r — (109) 

i=l  ^  ik=l 

with  the  coefficients  Ck  defined  by  the  formula 


ck  =  wkYlqie-^'‘^"~"^\ 

i=l 


(110) 
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for  all  k  =  1,2,. . .  ,No.  We  will  view  (109)  as  an  expansion  of  the  function  $  into  a  linear 
combination  of  exponentials  centered  at  w.  Now,  given  another  point  w  e  Bi, 

(109)  assumes  the  form 

No 

~  (111) 

fc=i 


and,  obviously, 

=  ),  (112) 
for  eill  k  =  1,2,. .. .  Thus,  we  are  lead  to  the  following  observation. 


Observation  3.20  An  expression  of  the  form  (111)  represents  the  potential  of  an  ar¬ 
bitrary  combination  of  charges  located  in  the  region  and  the  representation  is  valid 
for  all  z  E  Bx.  Furthermore,  (112)  can  be  interpreted  to  mean  that  in  the  representa¬ 
tion  (111),  the  translation  operator  for  the  potential  fields  is  diagonal.  In  the  following 
subsection  we  formalize  this  observation. 


3.5  Detailed  Description 

In  this  section,  we  use  Lemma  3.19  to  obtain  exponential  representations  of  the  form 
(111)  for  potentials  generated  by  collections  of  charges  located  in  Q*"  at  all  scales  s  >  0. 
Error  estimates  for  such  representations  are  provided  by  Theorem  3.21  below.  We  start 
with  the  following  obvious  generalization  of  Lemma  3.19. 


Theorem  3.21  Suppose  that  e  >  0  is  a  real  number,  and  real  numbers  xi,X2, . . .  ,xn^ 
and  wi,W2, . . .  fWNo  are  such  that 


No 

Wn 

n=l 


<  e 


(113) 


for  all  zq  e  0‘”,  z  e  Bi.  Suppose  further  that  s  >  0  is  a  real  number,  w  G  Bs,  and  the 
functions  En,s,w  n=l, 2,. ..  ,Nq ),  are  defined  by  the  formula 


En,s,w{z)  =  -  exp(-x„ 
s 


for  all  z  E  C.  Finally,  suppose  that 


(114) 


m 

=  (115) 

is  the  potential  due  to  a  set  of  m  charges  of  strengths  9i,?2, located  at  points 
zi,Z2,...,Zm  inside  the  square  L!'".  Then  for  any  z  e  Bs 

No  .  m 

l^(«)  -  £  bnEn,s,Uz)\  <  -  £  kil,  (116) 

n=l  ^  i=l 
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K  =  exp(-i„— - —), 


(117) 


for  all  n  =  1,2,. . .  ,Nq. 


The  following  theorem  is  an  immediate  consequence  of  the  definition  (114)  of  the 
functions  En,s,w 

Theorem  3.22  Suppose  that  under  the  conditions  of  Theorem  3.21,  w,w  E  Bs  and 


$(2)  ~  y 


(118) 


n=l 

are  the  expansions  provided  by  Theorem  3.21.  Then 


k  (  w-w 

bn  =  eXp{-Xn - )  Cr 


(119) 


(120) 


for  all  n  =  1,2,. . .  ,No. 


Suppose  that  s  >  0  is  a  real  number,  w  E  Bg  and  En,s,w  are  functions  defined 

by  (58)  and  (114),  respectively.  We  will  define  coefficients  rn,k,s,w  hy  the  formula 


^  t  _ _ 

^n,k,s,v}  =  ~  /  ^fc)^s,o(0  ^n,s4{w)  ds{^), 


(121) 


where  si,  S2,...,  are  given  by  (37)  and  k  =  1,2,. . . ,  and  n  =  1,2,. . . ,  Nq. 


Theorem  3.23  Suppose  that  s  >  0  is  a  real  number,  $  given  by  the  formula  (63)  is  the 
potential  due  to  a  set  of  m  charges  of  strengths  qi,q2,. . .  ,qm  located  in  the  square  Q®", 
and 


4(2)  =  £  a„«“*„(2), 


(122) 


is  the  expansion  provided  by  Theorem  3.8.  Suppose  further,  that  e  >  0  is  a  real  number 
and 

No 

$(2r)  -  bnEn,sA^),  (123) 

n=l 

is  the  approximation  provided  by  Theorem  3.21.  Then  for  n  =  1,2,.  ..,No 


^  ^  f'n,k,s,‘W  (^k^ 


(124) 


with  rn,k,s,w  defined  by  (121). 
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Proof.  Due  to  (58),  Theorem  3.1  and  Theorem  3.21 


nioM  =  «r‘(f )  =  «r(f ) 

1  iVo  ,  _ 

~  7=  5^'  E 

V*  n=l 

=  7?  E  <»««“'"'“>  /  5-'e-*-(=T‘i  fr(^)*(f) 

V<S  „=1  -/jOJ"  5 


No 


f'n,k,s,wEn,s,wiz)t 


71=1 


(125) 


for  each  A:  =  1,2, - Now  (124)  follows  immediately  from  (125). 


□ 


Suppose  that  s  >  0  is  a  real  number,  w  £  and  En,s,w  are  functions  defined 

by  (58)  and  (114),  respectively.  We  will  define  coefficients  en,k,s,w  by  the  formula 


^n,k,s,w  —  En,s,w{C} 


(126) 


k  =  1,2,. . . ,  and  n  =  1,2,. . . ,  Nq. 

The  proof  of  the  following  theorem  is  quite  similar  to  that  of  Theorem  3.8,  and  is 
omitted. 


Theorem  3.24  Suppose  that  s  >  0  is  a  real  number,  $  given  by  the  formula  (63)  is  the 
potential  due  to  a  set  of  m  charges  of  strengths  9i,925---j9m  located  in  the  square 
and 

OO 

»(^)  =  E  (127) 

Jk=l 

is  the  expansion  provided  by  Theorem  3.9.  Suppose  further,  that  e  >  0  is  a  real  number 
and 

No 

$(z)  ~  anEn,sA^),  (128) 

n=l 

is  the  approximation  provided  by  Theorem  3.21.  Then  fork  =  1,2,. . . 

No 

^k  ^  ^  ^n,k,s,w  On,  (129) 

n=l 


“with  Cn^kySyW  defned  by  (126). 


23 


Level  0 


Level  1 


Level  2 


Level  3 


Figure  3:  The  computational  box  and  three 
levels  of  refinement. 


4  The  multipole  algorithm 

4.1  Notation 

Without  a  loss  of  generality  we  can  assume  that  all  particles  are  located  in  the  unit 
square  centered  at  the  origin.  We  will  refer  to  this  square  as  the  computational  box.  A 
hierarchy  of  meshes  is  introduced  in  the  computational  box.  Mesh  level  0,  denoted  by 
5o,  corresponds  to  the  entire  computational  box.  Mesh  level  I  +  1,  denoted  by  5/+i, 
is  obtained  from  mesh  level  I  by  subdividing  some  boxes  into  four  equal  squares  (see 
Figure  3),  which  will  be  referred  to  as  children  of  the  given  square.  We  fix  an  integer 
number  s  >  0  and  at  each  level  we  subdivide  only  those  boxes,  which  contain  more  than 
s  particles.  A  box  which  is  not  subdivided  is  called  childless. 

Colleagues  of  a  box  are  adjacent  boxes  of  the  same  size. 

Two  boxes  of  the  same  size  which  are  not  adjacent,  are  called  separated. 

For  each  box  b  there  are  four  lists  of  other  boxes,  defined  as  follows. 

List  1  of  a  box  b  is  denoted  by  t4;  if  b  is  childless,  it  consists  of  b  and  all  childless 
boxes  adjacent  to  b.  If  6  is  a  parent  box,  then  £4  is  empty. 

List  2  of  a  box  b  is  denoted  by  H.  Vb  is  formed  by  all  the  children  of  the  colleagues 
of  fe’s  parent  which  are  separated  from  b  (see  Figure  4). 

List  3  of  a  box  b  is  denoted  by  W^,.  If  b  is  childless,  VFj,  consists  of  all  descendants  of 
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Figure  4:  Interaction  List  2  for  a  dashed  box. 


6’s  colleagues  whose  parents  are  adjacent  to  6,  but  which  are  not  adjacent  to  h  themselves. 
If  6  is  a  parent  box,  Wb  is  empty. 

List  4  of  a  box  b  is  denoted  by  Xy,  it  consists  of  all  boxes  b'  such  that  6'  €  Wj,. 

will  denote  the  p-term  expansion  of  the  form  (64)  of  the  potential  due  to  particles 
located  in  b. 

will  denote  the  p-term  expansion  of  the  form  (64)  of  the  potential  due  to  particles 
located  outside  UbUWt. 


4.2  Informal  description  of  the  algorithm 

The  data  structure  used  by  our  algorithm  is  virtually  identical  to  the  one  presented  in 
[2].  It  relies  on  clustering  of  particles  at  various  scales.  The  computation  of  interactions 
between  clusters  separated  from  each  other  is  performed  via  the  expansions  (64)  and 
(67),  while  the  interactions  between  nearby  particles  are  computed  directly. 

In  order  to  adapt  our  grid  to  a  given  distribution  of  particles  we  fix  an  integer  s  >  0 
and  subdivide  only  those  boxes  which  contain  more  than  s  particles. 

The  algorithm  consists  of  the  following  stages. 

(1)  We  create  a  hierarchy  of  meshes  in  a  computational  cell. 

(2)  For  each  childless  box  b  we  directly  evaluate  interactions  between  particles  in  b 
and  particles  in  Ub,  List  1  of  6. 

(3)  For  each  childless  box  b  we  form  the  expansion  into  outgoing  singular  func¬ 
tions  by  means  of  Theorem  3.8. 
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(4)  For  each  parent  box  b  we  form  the  expansion  by  merging  expansions  of  its 
children  via  Theorem  3.12. 

(5)  We  convert  all  expansions into  the  exponential  form  via  Theorem  3.23. 

(6)  We  use  Theorem  3.22  to  shift  the  q-term  exponential  expansion  to  each  box  in 
Vb,  List  2  of  b. 

(7)  For  every  particle  r  €  6,  we  compute  the  field  due  to  all  particles  in  Wt,  List  3  of 
b,  by  summing  up  the  expansions  for  all  w  eWb  and  add  it  to  the  potential  at  this 
point  calculated  in  (1). 

(8)  We  convert  the  field  of  each  particle  in  Xb,  List  4  of  6,  into  the  expansion 

(9)  We  convert  all  exponential  expansions  into  expansions  via  Theorem  3.24  and 
combine  them  with  the  result  of  (8). 

(10)  For  each  child  box  b  we  shift  the  expansion  of  its  parent  (Theorem  3.9)  and  add 
it  to  the  expansion  ^5”. 

(11)  For  each  childless  box  b  we  evaluate  the  expansion  at  every  particle  r  €  6 
and  add  it  to  the  result  of  (7),  obtaining  the  field  at  r. 

4.3  Detailed  description  of  the  algorithm 

ALGORITHM 

Comment  [Choose  main  parameters.] 

Choose  precision  e  to  be  achieved.  Set  the  length  of  expansions  according  to 
Table  1. 

Set  the  maximum  number  s  of  particles  in  a  childless  box. 

Step  1 

Comment  [Refine  a  computational  cell  into  a  hierarchy  of  meshes.] 

do  1  =  1,2,. . . 

do  bi  €  Bi 

if  bi  contains  more  than  s  particles  then 

subdivide  bi  into  four  boxes  and  add  the  nonempty  boxes 
formed  to  Bi+i. 
end  if 
end  do 
end  do 

Comment  [We  denote  by  nlev  the  highest  level  of  refinement  and  by  nbox  the  total 
number  of  boxes  formed  at  Step  1.] 
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step  2 


Comment  [For  every  particle  in  each  childless  box  h  compute  directly  the 
interactions  with  particles  in  t/j,,  List  1  of  6.] 

do  i  =  l,nbox 

if  hi  is  a  childless  box  then 

for  each  particle  r  Eh  compute  interactions  between  r  and 
all  particles  in  Ub^. 
end  if 
end  do 


Step  3 

Comment  [For  every  childless  box  b  form  an  expansion  into  outgoing  singular 
functions  about  the  center  of  b  via  Theorem  3.8.] 

do  i  =  l,nbox 

if  bi  is  a  childless  box  then 

use  Theorem  3.7  to  form  a  p-term  expansion 
representing  the  potential  due  to  all  charges  in  h. 

end  if 
end  do 


Step  4 

Comment  [For  each  parent  box  b  use  Theorem  3.10  to  shift  the  center  of  each  6’s 
child  box’s  expansion  to  6’s  center  and  add  the  resulting  expansions 
together.] 

do  1  =  nlev-1,1,-1 
do  bi  €  Bi 

if  bi  is  a  parent  box  then 

use  Theorem  3.10  to  obtain  a  p-term  expansion  by 
shifting  expansions  of  Bi's  children  to  b  and  adding  the 
resulting  expansions  together. 

end  if 
end  do 
end  do 


Step  5 
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Comment  [For  each  box  b  use  Theorem  3.23  to  convert  the  expansion  to  the 
exponential  form.] 

do  i  =  l,nbox 

use  Theorem  3.23  to  convert  the  expansion  to  the  exponential  form, 
end  do 


Step  6 

Comment  [For  each  box  h  use  Theorem  3.22  to  shift  the  exponential  expansion  from 
b  to  each  box  in  Vb,  List  2  of  6.] 

do  i  =  l,nbox 

use  Theorem  3.22  to  shift  the  exponential  expansion  from  b  to  each  box 
in  Vb,  List  2  of  b. 
end  do 


Step  7 

Comment  [For  each  childless  box  b,  evaluate  the  expansions  for  all  w  €  Wb, 
List  3  of  b,  at  every  particle  r  €  b.] 

do  i  =  l,nbox 

if  b  is  childless  then 

evaluate  the  expansion  of  each  box  tn  €  VFj  at  every  particle 
r  £  b. 

end  if 
end  do 


Step  8 

Comment  [For  each  box  b,  create  the  expansion  'S'j”  representing  in  b  the  field  due  to 
particles  in  Xb,  List  4  of  b.] 

do  i  =  l,nbox 

Create  the  expansion  of  the  field  due  to  all  particles  in  Xb. 
end  do 


Step  9 

Comment  [For  each  box  b  use  Theorem  3.24  to  convert  the  exponential  expansion 
into  the  expansion  of  the  form  (67)  and  combine  with  the  expansion 
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computed  in  Step  7.] 
do  i  =  l,nbox 

use  Theorem  3.24  to  convert  the  exponential  expansion  the  expansion  of 
the  form  (67)  and  computed  it  with  the  expansion 
end  do 


Step  10 

Comment  [Use  Theorem  3.12  to  shift  the  local  expansions  of  parent  boxes  to 
their  children.] 

do  1  =  l,nlev-l 
do  hi  €  Bi 

if  hi  is  a  child  box  then 

use  Theorem  3.12  to  shift  a  p-term  expansion 
from  6’s  parent  to  b  and  update 
end  if 
end  do 
end  do 


Step  11 

Comment  [For  every  childless  box  b  evaluate  the  expansion  '4'^”  at  each  particle 
r  E  b  and  add  it  to  the  result  of  Step  7,  obtaining  the  potential  at  r.] 

do  i  =  l,nbox 

if  bi  is  childless  then 

for  each  particle  r  E  bi  evaluate  the  p-term  expansion 
obtaining  the  potential  at  location  r. 
end  if 
end  do 

4.4  Complexity  analysis 

Step  Operation 
count 

1  Np  Each  particle  is  assigned  to  a  box  at  every  level.  There  are 

at  most  p  levels  of  refinement. 

2  33Nps  The  direct  computation  of  interactions  between  any  two 
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3 

4 

5 

6 
7 


Np 


SNp^ 

s 


SNp^q 

s 

135Npq 

$ 


32iVp2 


childless  boxes  requires  at  most  ^  operations.  The  total 
number  of  boxes  appearing  on  a  List  1  of  a  box  does  not 
exceed  (see  [2]). 

Each  paxticle  contributes  to  a  p-term  expansion  of 
a  unique  childless  box  b  containing  it. 

Each  translation  requires  operations  and  there  are  at  most 
boxes  (see  [2]). 

Each  translation  requires  pq  operations.  The  total  number 
of  boxes  is  bounded  by 

Each  diagonal  translation  requires  q  operations  and  there 
are  at  most  27  boxes  on  any  List  2. 

Computing  the  interactions  of  all  particles  in  a  box  b  with 
a  box  in  Wb  requires  ps  operations.  The  total  number  of 
boxes  appearing  on  a  List  3  does  not  exceed  (see  [2]). 


8  Z2Np^  Computing  the  interactions  of  a  box  b  with  all  particles  in 
a  box  in  Xb  requires  ps  operations.  The  total  number  of 
boxes  appearing  on  a  List  4  does  not  exceed  (see  [2]). 


9 

10 
11 


SNp^q 

s 


SNp^ 

$ 


Np+N 


Each  translation  requires  pq  operations.  The  total  number 
of  boxes  is  bounded  by 
Each  translation  requires  p^  operations. 

A  p-term  expansion  is  evaluated  at  each  particle. 
Summing  up  requires  extra  N  operations. 


Combining  the  CPU  times  for  all  the  above  stages,  we  obtain  the  estimate 

3  2 

T  =  N{a  —  +  b^  +  cps  +  dp^),  (130) 

s  s 

where  the  constants  a,  6,  c,  d  depend  on  the  implementation.  The  parameter  s,  the  maxi¬ 
mum  number  of  particles  in  a  childless  box,  should  be  chosen  so  as  to  minimize  the  CPU 
time.  An  elementary  calculation  shows  that  the  minimal  time  Tmin 


is  obtained  for 


Tmin  =  ATp^^ 

1(^  +  0-, 

1  p 

(131) 

5  • 

^mtn  —  f-\ 

y/c\ 

p 

(132) 
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with  constants  a,  ^  dependent  on  the  particular  implementation.  Since  p  ~  we  arrive 
at  the  estimate 

Tmin  =  'yNp^  =  -yN  (log2  (133) 

where  the  constant  7  depends  on  the  implementation. 

The  storage  requirements  are  determined  by  the  lentgh  of  expansions  and  the  total 
number  of  boxes,  which  does  not  exceed  Per  box,  we  store  2p  coefficients  of  singular 
functions  and  4^  coefficients  of  exponentials.  Additionally,  we  store  the  position  and  the 
charge  of  each  particle  and  corresponding  p  values  of  singular  functions.  Therefore,  the 
storage  requirements  are  of  the  form 


5  =  iv(c,  + 

s  ^ 

with  the  constants  ci,C2  dependent  on  the  implementation.  Ks  =  Smin,  we  obtain 

5  =  C3  iVlogj  J, 

where  C3  depends  on  the  implementation. 


(134) 


(135) 


Remark  4.1  A  careful  examination  of  (135)  shows  that  even  though  the  storage  re¬ 
quirements  of  the  algorithm  are  proportional  to  the  number  of  particles,  the  associated 
constant  is  quite  large,  especially  in  double  precision  calclations.  This  limits  the  size  of 
problems  which  can  be  handled  in  computing  environments  where  the  available  memory 
is  limited.  Moreover  in  systems  with  virtual  memory,  it  is  liable  to  increase  the  wall-clock 
time  of  the  algorithm.  This  problem  is  presently  being  addressed  by  the  authors. 


5  Numerical  results 

The  algorithm  described  in  Section  4  has  been  implemented  in  Fortran  77  and  numerical 
experiments  have  been  performed  on  a  SPARCstation  2.  We  compare  its  performance 
with  an  implementation  of  the  FMM  from  [2],  and  with  direct  application  of  the  poten¬ 
tial  matrix.  We  give  the  results  for  three  regimes:  particles  uniformly  distributed  in  a 
square,  particles  located  uniformly  on  a  curve  and  particles  clustered  within  a  square. 
All  calculations  were  done  in  double  precision,  and  extended  (quadruple)  precision  was 
used  to  determine  the  relative  errors.  The  number  of  particles  varied  between  400  and 
25,600,  with  charge  strengths  randomly  distributed  on  the  interval  (0, 1). 

The  results  of  our  experiments  are  given  in  Tables  (6-8).  In  each  table  the  first 
column  contains  the  number  N  of  particles  in  a  simulation.  Second, third  and  forth 
columns  show  the  CPU  times  Tnew^Toid^Tur  in  seconds  of  the  present  algorithm,  the 
algorithm  described  in  [2]  and  of  the  direct  calculation,  respectively.  The  fifth,  the  sixth 
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Table  6:  Uniformly  distributed  particles. 


N 

T 

new 

Told 

Tdir 

Enew 

Eoid 

Edit 

400 

0.3 

1.5 

1.0 

2.2  10-1® 

4.4 10-1® 

4.3 10-1® 

800 

0.8 

2.7 

4.1 

3.410-1® 

1.4  10-1® 

6.110-1® 

1600 

1.5 

7.1 

16.9 

3.4 10-1® 

5.2  10-1® 

2.2  10-1® 

3200 

3.5 

10.9 

68.6 

3.110-1® 

1.3  10-1® 

8.5 10-1® 

6400 

6.4 

32.2 

277.3 

2.410-1® 

5.6  10-1® 

1.310-1® 

12800 

12.9 

46.9 

(1100) 

3.610-1® 

4.9  10-1® 

— 

25600 

26.7 

143.5 

(4400) 

4.9 10-1® 

4.6  10-1® 

— 

Table  7:  Particles  distributed  on  a  curve. 


N 

T 

new 

Told 

Tdir 

Enew 

Eoid 

Edir 

400 

0.3 

1.2 

1.0 

2.6 10-1^ 

1.7  10-1^ 

2.5  10-1^ 

800 

0.6 

2.5 

4.2 

2.3 10-1^ 

3.210-1'! 

2.3  lO-i'i 

1600 

1.2 

3.8 

17.2 

7.7 10-1^ 

9.2  10-1^ 

5.110-1^ 

3200 

2.6 

7.7 

69.4 

3.4 10-1^ 

4.1 10-1^ 

8.9 10-1^ 

6400 

5.4 

15.4 

283.2 

2.1 10-1® 

2.8  10-1® 

3.110-1® 

12800 

10.9 

30.7 

(1100) 

2.3 10-1® 

3.7  10-1® 

25600 

21.6 

58.8 

(4400) 

7.4 10-1® 

6.9  10-1® 

— 

and  the  seventh  column  contain  the  corresponding  relative  errors  EnewiEoid^Enr.  The 
errors  are  computed  in  the  P  norm  via  the  formula 


(Ef=. 


(136) 


where  is  the  value  of  the  potential  at  the  ith  particle  position  obtained  in  a  direct 
calculation  in  extended  (quadruple)  precision,  and  $2*  is  the  corresponding  value  obtained 
by  one  of  the  three  methods  in  question. 

In  the  first  set  of  experiments  particle  positions  are  generated  randomly,  resulting 
in  their  nearly  uniform  distribution  (Figure  5).  In  the  second  set  of  tests  particles  are 
uniformly  distributed  on  an  ellips  (Figure  6).  The  third  set  of  tests  is  performed  for  a 
nonuniform  distribution  of  particles  within  a  square  (Figure  7). 

In  all  cases  the  maxium  number  of  particles  in  a  childless  box  is  set  to  44. 


33 


Figure  7:  Particles  distributed  on  a  curve. 
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Table  8:  A  nonuniform  distribution  of  particles. 


N 

T 

J-new 

Told 

Tdir 

Enew 

Eoid 

Edir 

400 

0.5 

1.4 

1.0 

3.5  10-^3 

3.5  10-12 

3.0 10-12 

800 

1.2 

2.8 

4.2 

1.110-^2 

1.110-12 

1.110-12 

1600 

2.3 

5.2 

17.1 

1.110-^2 

1.1 10-12 

1.110-12 

3200 

4.9 

11.6 

69.0 

3.910-^3 

3.9  10-12 

3.710-12 

6400 

8.7 

25.8 

285.2 

4.9  10-12 

4.9  10-12 

2.710-12 

12800 

17.8 

52.6 

(1100) 

6.710-12 

6.8  10-12 

— 

25600 

34.7 

110.6 

(4400) 

2.510-12 

2.510-12 

— 

The  following  observations  can  be  made  from  the  experiments: 

(1)  The  CPU  time  of  the  present  algorithm  grows  linearly  with  N.  The  breakeven 
point  with  the  direct  calculation  ranges  from  No  =  140  for  uniformly  distributed  particles 
to  No  =  190  for  a  highly  nonuniform  distribution. 

(2)  The  performance  of  the  algorithm  does  not  depend  significantly  on  the  type  of  the 
particle  distribution.  The  timings  for  a  highly  nonuniform  distribution  are  about  60% 
higher  that  in  the  case  of  uniformly  distributed  particles. 

(3)  The  accuracy  obtained  by  the  algorithm  agrees  with  the  error  bounds  (3.15)  and 
(3.13). 

6  Conclusions 

A  new  version  of  the  Fast  Multipole  algorithm  for  Coulombic  interactions  has  been  devel¬ 
oped.  While  the  prior  schemes  use  Laurent  and  Taylor  series  to  represent  the  potential, 
we  use  singular  functions  of  an  appropriately  chosen  operator,  obtaining  a  much  faster 
convergence.  We  also  introduce  an  intermediate  representation,  in  which  most  transla¬ 
tion  operators  are  diagonal.  In  two  dimensions,  the  resulting  scheme  is  three  to  five  times 
faster  than  the  best  implementations  of  the  old  FMM  we  are  aware  of;  we  expect  the 
ratio  to  be  much  higher  in  three  dimensions. 

The  CPU  time  estimate  and  the  storage  requirements  for  the  algorithm  are  of  the 
order  of  0{N),  where  N  is  the  number  of  particles.  Its  performance  does  not  depend 
significantly  on  the  particle  distribution. 

The  three-dimensional  version  of  the  algorithm  is  being  implemented,  and  will  be 
presented  at  a  later  date. 
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Appendix  A 


Estimates  for  the  singular  values  of  the  operators  Cm 
and  Cout 


In  this  section  we  will  express  the  operators  ^ 

Cout  '■  L^{D^’^)  by  means  of  the  Cauchy  integral  operator  commonly  encoun¬ 

tered  in  complex  analysis  and  derive  inequalities  (39)  from  corresponding  inequalities  for 
the  new  operator. 

We  will  define  arclength  parametrizations  71  ;  [0,4]  Z)*”,  72  :  [0, 12]  of  the 

squares  respectively,  by  the  formulas 


7i(<)  = 


\-t  +  ¥ 

for 

0  <t  <  1 

for 

1  <t  <  2 

-k  +  t-li 

for 

2  <  t  <  3 

2  +  (“2  + 

for 

3  <  t  <  4 

(137) 


72(i)=3  7i(|)  for  0<t<12.  (138) 

Obviously,  except  for  the  corners  of  the  squares,  71,  72  are  differentiable  and 

KI  =  KI  =  i. 

Finally,  we  define  the  Cauchy  integral  operator  Cq  :  ->  by  the 

formula 


Co/(^)  = 


1 

27r? 

1 

27ri 


/(O 


( m 

hi- 


di 


ti-z 

mm 


L 


Din  ^  —  Z 


MO- 


(139) 


Observation  A.l  Denoting  hy  unitary  operator  of  mul¬ 

tiplication  by  Yi,  we  observe  that  Co  =  CoutM^[-  Therefore,  the  singular  values  of  the 
operator  Cout  o,i'e  identical  to  those  of  Cq. 


Suppose  that  ri,  r2  are  real  numbers  such  that  ^ 


<  ri  <  r2  < 


|.  We  will 


denote  by  Si  and  S2  the  circles  =  {|2r|  =  rj}  and  ^2  =  {l^j  =  r2}  oriented  coun 
terclockwise.  We  define,  respectively,  the  integral  operators  Cp  :  L\S2)  Z2(Z>°“‘), 

:  L\Si)  ->  L\S2),  CY  :  Z2(Z>-)  L\Si)  by  the  formulae 


(141) 


^2 


f(0 


d(, 


(142) 


Lemma  A. 2  If  the  operators  Cq,  Cp,  Cp  are  defined  by  formulae  (139-142),  re¬ 

spectively,  and  ^  <  Ti  <  r2  <  |,  then 

Co  =  (143) 


Proof.  Applying  twice  Cauchy’s  integral  formula  for  an  unbounded  domain  (see  [6,  vol.l, 


p.318]),  for  an  arbitrary  function  /  €  L^(£>’”)  we  have 

C?’^^{C;^f){0  =  iC  €  52),  (144) 

C?{C?f){z)  =  -c;^f{z),  (z  e  D°“‘).  (145) 

Combining  (144)  with  (145),  we  obtain 

Cp(Cp’^^Cp/)(z)  =  Cp(-Cp/)(^)  =  Cp/(^)  =  Cofiz),  (146) 

for  any  z  €  D°^K  □ 


The  following  lemma  is  well-known  (see  [3,  p.98  and  p.l44]). 

Lemma  A. 3  Suppose  that  Xi,  X2,  X3,  X4,  are  piecewise  smooth  curves  in  R,^.  If 
A  :  L^iXz)  L^(Xi),  B  :  L^iX^)  1^X2),  C  :  1^X4)  L^X^)  are  bounded  in¬ 

tegral  operators,  and  B  is  a  Hilbert-Schmidt  operator,  then  the  composition  ABC  is  a 
Hilbert-Schmidt  operator.  Moreover,  if  we  denote  by  Sn{B)  and  Sn{ABC)  the  singular 
values  of  the  operators  B  and  ABC,  respectively,  then 

s„(A5C)<||A||  lie'll  s„(.B),  (147) 


for  all  n=l,2, _ 

Lemma  A. 4  Suppose,  that  {^n}  and  {V’n}  o,re  orthonormal  systems  in  L^{Si)  and 
L^{S2),  respectively  defined  for  n=l,2,. by  formulae 


MO  = 

y/2-Kri  ( 

!)■ 

(148) 

MO  =  ■ 

1 

(149) 

\/27rr2 

\z)  • 
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Then 


= f:  (-)"'*  ( /,  ^»>  v>. .  (150) 

n=l  ^’"2/ 

IS  a  Singular  Value  Decomposition  for  the  operator  :  L^{Si)  —*■  L'^(S2).  In  particu¬ 
lar,  the  singular  values  SniC^^'^^)  of  €2^'^^  are  given  by  the  formulae 

=  (^)"‘’  ,  (151) 

for  n=l,2,.... 


Proof.  Obviously,  for  any  ^  £  S\,  z  £  S2, 

1 


(152) 


Combining  (141)  with  (152),  we  obtain 

-  1  f  m 1 


1  00  . 

1  i*  _ 

=  /  /(0\^’-r«f)*(0 

n=i  JSi 

=  £  (~)  ( /.  <^0  V’n(^)  • 

n=1  ^^2/ 


(153) 


The  following  two  theorems  establish  bounds  from  above  and  from  below  for  the 
singular  values  of  the  operator  Cout-  We  only  prove  Theorem  A.5  below;  the  proof  of 
Theorem  A.6  is  nearly  identical,  and  is  omitted. 

Theorem  A.5  For  any  real  number  ^  <  p<l  there  exists  a  constant  c>  0,  such  that 

^n(fiout')  ^  Cp  j  (154) 

for  all  n=l,2,...,  where  Sn{Cout)  are  the  singular  values  of  the  operator  Cout- 
Proof.  Defining  real  numbers  ri,  r2  by  the  formulae 


ri 


\/2  +  3/3 
4 


(155) 


40 


(156) 


^1 

r2  =  — 

P 

we  observe  that 

y/2  S  ,  , 

—  <  <  ^2  <  2’  (157) 

Let  us  consider  the  operators  Co,  Cp,  Cl^  defined  by  formulae  (139-142).  Clearly, 

^  <  Ti  <  r2  <  |.  p=  Combining  (143)  with  (147),  we  obtain 

sn{Co)  =  <  \\cn  ||.  (158) 

According  to  Observation  A.l,  the  singular  values  of  operators  Cout  and  Cq  coincide.  Now 
Lemma  A.4  gives 

<s.(C<)  =  s.(Co)  <  licni  Sn(c;‘''’)  IICJMI  =  c  (n)“  =  ep",  (159) 

for  all  n=l,2,. . . .  □ 


Theorem  A.6  For  any  number  0  <  p  <  there  exists  a  constant  c  >  0,  such  that 

Sn{CoutT>  cp^,  (160) 

for  all  n=l,2,...,  where  Sn{Cout)  denote  the  singular  values  of  the  operator  Cout- 
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Appendix  B 


Estimates  for  the  singular  functions  of  the  operators 

^in  und  Cout 


In  this  section  we  estimate  the  maxima  of  singular  functions  of  the  operator  Cout-  We 
begin  with  four  technical  lemmas. 

Lemma  B.l  Ifn  >  0  is  an  integer  and  P  :  €  ^  €  is  a  polynomial  of  degree  n-1,  then 

^  ||-P||l2[-i,i]  5  (161) 

(for  definitions  of\\  •  ||l2[-i,i]  and  ||  •  ||l®°[-i,i]  see  Section  2.1). 

Proof.  Let  Pk  denote  the  Legendre  polynomials  on  the  segment  [-1, 1].  It  is  well-known 
(see  [1]),  that 

ll-ffcllL^hi,!]  =  1,  (162) 

(163) 

Obviously,  there  exist  unique  complex  coefficients  co,ci, . . .  ,Cn_i  such  that 

n-1 

■P  =  Cfc  Pk,  (164) 

k=0 

and,  consequently. 


n— 1 


n— 1  c\ 

2  ^ 


=  E  lc»P  llftifei-.,.)  =  E  I'^il 

k=o 

Now,  due  to  Schwarz’s  inequality,  we  have 


k=o  jb=o  2A:  -1- 1 


(165) 


n— 1 


n— 1 


||P|U“[-1,1]  <  Yj  l^fel  l|Pfe|U°°[-l,l]  =  S  l^fel  =  ^  ^  kfcl 

A:=0  A:=0  k=0  V  ^  V  2.K:  +  1 


^n-1 


2  E 


2  /n— 1 


\k=0 


n 


||P||l2[-i,i]-  (166) 

□ 


The  following  lemma  is  readily  obtained  from  Lemma  B.l  by  scaling. 
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Lemma  B.2  Suppose,  that  n  >  0  is  an  integer  and  P  :  C  — >  C  is  a  polynomial  of  degree 
n  —  1.  If  I  C  ^  is  any  segment  of  length  a  >  0,  then 


||-P||l«(7)  <  \\P\\l2(I)-  (167) 

Lemma  B.3  There  exists  a  constant  Ci  >  0  such  that  for  any  segment  I  C  of  length 
a  >  0,  integer  n  >  0  and  any  f  €  there  exists  a  polynomial  P  of  degree  n  —  1 

satisfying 

W^outf  —  <  Cl  o”  II/II25  (168) 

■with  the  operator  Cout  •  L'^[D°^^)  defined  by  (30). 

Proof.  DiiFerentiating  (30)  n  times  we  obtain 

(169) 

and,  due  to  Schwarz’s  inequality,  we  have 

ll(C<../)<">|U<c,n!  ll/lb,  (170) 

where  Ci  >  0  is  a  constant.  Let  us  denote  by  zq  the  midpoint  of  the  segment  I  and  let  P 
be  the  Taylor  polynomial  of  order  n  —  1  centered  at  zq  for  the  function  Coutf-  Now,  the 
Taylor’s  formula  implies  that 

|C~./(^)  -  />(^)l  <  I(c„,/)<’>(C)I,  (171) 

and  therefore 

l|C„,/  -  PIIl-i,)  <  ^  ||(C„,/)(">|U  <  Cl  ^  ll/lb  <  Cl  a"  ll/lb.  (172) 

□ 


Lemma  B.4  There  exists  a  constant  c  >  0;  such  that  for  any  f  G 

l|C..,/|U  <  c ||C„,/|b  (log +  l)  . 

V  WtoutJh  J 

where  the  operator  Cout  •  L^{D''^)  is  defined  by  (30). 

43 


(173) 


Proof.  Let  0  <  a  <  1  be  a  real  number  and  I  any  segment  of  length  a  included  in  the 
square  For  an  integer  n  >  0  and  /  €  P{D^)  denote  by  P  the  polynomial  provided 
by  Lemma  A. 9.  Due  to  Lemma  A.8  we  have 


<  \\Coutf  —  ^  ||-P|U2(/). 


(174) 


We  will  need  the  following  obvious  inequality 


<  \\Coutf  —  P||l*(/)  +  \\Coutf\\L^(I) 

<  Vg  \\Coutf  —  P||l<»(/)  +  \\Coutf\\L^(I)' 


(175) 


Combining  (168),  (174)  and  (175)  we  obtain 


\Koutf\\Lo^{I)  <  (»^  +  1)  \\Coutf  —  PIlL^g)  +  ^  ||<^out/||L2(/) 


<  C2  n  {a^Wfh  +  WCoutfhj 

V  I|C~!/I|2  Va/ 


(176) 


with  some  constants  C2  >  0,  C3  >  0.  Setting  o  =  ■^  (where  e  is  the  base  of  natural 
logarithms),  and 


L  \\^outj\\2  _ 


we  arrive  at 


I|Co../IIl-(z)  <  e  ||C..,/||2  log  +  1 


(177) 


(178) 


where  c  >  0  is  a  constant.  Since  every  point  in  is  contained  in  some  segment 
I  Q  £)out  q£  length  CL  =  \i  the  lemma  follows.  □ 

Now  we  proceed  to  the  principal  result  of  this  section. 

Theorem  B.5  Let  be  the  left  singular  functions  of  the  operators  Cin  andCout, 

respectively.  Then  there  exist  constants  0  <  q  <l  and  c>  0  such  that 


IICIIco  <  cn, 

UrWoo  <  cn, 


(179) 

(180) 


for  all  n=l,2,. . . . 
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Proof.  Since  Cout^T*  =  Lemma  A.  10  implies  that,  there  is  a  constant  Cj  >  0, 

such  that 


=  Cl  Sn  (log  —  +  log  ||C(,«t||  +  l)  . 

\  Sn  J 

Now,  due  to  Theorem  A. 6,  there  exists  a  constant  c  >  0,  such  that 


(181) 


<  cn. 


(182) 


for  n=l,2, - It  is  easy  to  see,  that  Lemma  A.  10  holds  with  Cout  replaced  by 

and  therefore  this  proof  extends  to  the  case  of  the  left  singular  functions  □ 
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